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ABSTRACT 

The HZETRN code was used in recent trade studies for renewed lunar exploration and currently used in 
engineering development of the next generation of space vehicles, habitats, and EVA equipment. A new version of 
the HZETRN code capable of simulating high charge and energy (HZE) ions, light-ions and neutrons with either 
laboratory or space boundary conditions with enhanced neutron and light-ion propagation is under development. 
Atomic and nuclear model requirements to support that development will be discussed. Such engineering design 
codes require establishing validation processes using laboratory ion beams and space flight measurements in 
realistic geometries. We discuss limitations of code validation due to the currently available data and recommend 
priorities for new data sets. 


INTRODUCTION 

Our early interest in transport code development paralleled the experimental studies of Walter 
Schimmerling et al. (1986) at the Lawrence Berkeley Laboratory and involved solving the Boltzmann transport 
equation for mono-energetic ion beams in the context of the continuous slowing down approximation (Wilson et al. 
1984). The main computational limitation was found to be the inadequacy of available nuclear data through 
comparing computational results with the ionization data for a broad beam of 20 Ne ions (Wilson et al. 1984, Shavers 
et al. 1990, 1993). Emphasis was soon overtaken by the need to establish the scope of the GCR protection problem 
and marching procedures were used to get first order estimates of shielding requirements (Wilson et al. 1991). 
Testing the new computational marching model against atmospheric air shower data again pointed to the 
inadequacy of the available nuclear data (Wilson and Badavi 1986) and development of a semi-empirical nuclear 
model followed, leading to reduced computational errors (Wilson et al. 1987a, 1987b). The next dozen years 
emphasized spaceflight validation of the marching solution and nuclear model improvements (Badhwar et al. 1996, 
Shinn et al. 1998, Cucinotta et al. 1998, Wilson et al. 1995, 2006a, 2007, Hugger et al. 2003, Nealy et al. 2007). 
Advanced solution methods of the Boltzmann equation continued to develop (Wilson et al. 1994a,b) but only 
slowly after NASA interest in deterministic transport code and nuclear model development shifted in favor of 
Monte Carlo methods (Armstrong and Colburn 2001, Pinsky et al. 2001). Recent renewed interest within NASA for 
deterministic code development is giving new emphasis to improved solution methods (Wilson et al. 2004, 2005a, 
Tweed et al. 2005, 2006) with a renewed focus on nuclear modeling activity (Cucinotta et al. 1994, 1998, 2007, 
Norbury et al. 2007). As a result, recent developments mainly utilize the older semi-empirical NUCFRG2 model 
(Wilson et al. 1995, 2005b, Walker et al. 2005) with energy downshifts and momentum dispersion (Tripathi et al. 
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1994). Comparison of NUCFRG2 with other models and experiments are given by Zeitlin et al. (1997), Cucmotta 
et al. (2007), and Walker et al. (2005). As a result for the need to design and test new systems for future 
exploration in the low Earth orbit environment, there has been renewed interest in flight validation using 
International Space Station as a measurements platform (Hugger et al. 2003, Wilson et al. 2005b, 2007, Nealy et al. 
2007). In the current report, we will review the current status of transport code development with emphasis on 
future needs. 

DETERMINISTIC CODE DEVELOPMENT 

The relevant transport equations are the linear Boltzmann equations derived on the basis of conservation 
principles (Wilson et al. 1991) for the flux density (j)j(x, Q , E) for particle type j as 

Q 0 V(/)j(x, Q, E)= Zfojk( 0,0', E, E ') Mx,0',E') dO' dE ' - oj(E) 4j(x,Q,E) (1) 

where oj(E) and ojk(Q,Q',E,E) are the shield media macroscopic cross sections. The ojk(Q f Q',E,E') represent 


all those processes by which type k particles moving in 
direction Q with energy E (including decay processes). 
particular product, and the appropriate cross sections for 
are functions of the particle fields and may be included 
once the particle fields are known. The total cross 
section oj (E) with the medium for each particle type is 

oj (E) = oj y at (E) + oj yC i (E) + oj yV (E) (2) 

where the first term refers to collision with atomic 
electrons, the second term is for elastic scattering on the 
nucleus, and the third term describes nuclear reactions 
where we have ignored the minor nuclear inelastic 
processes (excitation). The corresponding differential 
cross sections are similarly ordered. Many atomic 
collisions (~10 6 ) occur in a centimeter of ordinary 
matter, whereas ~10 3 nuclear coulomb elastic collisions 
occur per centimeter, while nuclear reactions are 
separated by a fraction to many centimeters depending 
on energy and particle type. Solution methods first use 
perturbations based on the ordering of the cross sections 
with atomic interactions as the first physical 
perturbation with special methods used for neutrons for 
which atomic cross sections are zero. 

We rewrite equation (1) in operator notation by 


defining a vector array field function as 

0 = [(f)j(x, Q, E) ] (3) 

the drift operator as 

D = [Q*V] (4) 

and the interaction operator as 

I = [Zfo jk (Q,a\E,E') dO dE' - oj(E)] (5) 

with the understanding that / has three parts associated 
with atomic, elastic, and reactive processes as given in 
equation (2). Equation (1) is then rewritten as 

[D-I at -I el ]*0 = l r -0 (6) 


where the first two physical perturbation terms are 
shown on the left-hand side and are represented by 


direction Q f with energy E produce a type j particle in 
Note that there may be several reactions that produce a 
equation (1) are the inclusive ones. Exclusive processes 



Fig. 1 Range of ions in aluminum. 
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Fig. 2 FWHM of 49.1 MeV protons. 


( 7 ) 


diagonal operators. The first order physical perturbation for atomic processes is solution of 

[D-IJ»9= 0 

using the moments methods, the solution can be approximated by 

9 (z,E) = exp[- (E - <E(z)>) 2 /(2 s 2 (z))]/(V(2x) s(z)) (8) 

where the array of mean residual energies <E(z)> and the energy deviation s(z) are evaluated using a second order 
Green’s function (Wilson et al. 2002) and related to range and full width at half maximum (FWHM) shown in Figs. 
1 and 2. Although straggling correction for the uncollided beam is important, it is negligible in the higher order 
terms compared to the energy dispersion in fragmentation. 

The second order physical perturbation is the Coulomb scattering by the atomic nucleus as a solution of 

[D-I el ]'*=0 (9) 

and represented by Rutherford scattering modified by 
screening of the nuclear charge by the orbital electrons using 
the Thomas-Fermi distribution for the atomic orbitals. We 
will utilize the multiple Coulomb scattering solutions of 
Fermi given by 

(z,r, 0 r ) = [V3 w 2 /2jvz 2 ] 

x exp[ - w 2 ( 0 2 /z-3r 6Jz 2 +3r 2 /z 3 ) ] (10) 

where z is the longitudinal distance, r is the lateral distance, 

8 r is the angle to the longitudinal axis, and w 2 the array of 
appropriate diffusion coefficients. Strictly speaking, the 
solution applies only over intervals for which the variation 
in ion energy is small. It follows that the mean square angle 
(understood as a differential quantity) in traveling a distance 
dz is given from equation (10) as 

( 0r) = 2dz/w 2 (11) 

Conversely, one finds for a uniform nuclear charge 
distribution shielded by a Thomas-Fermi atomic structure 

(dr) = dz (ZpEJfipcf/Xo (12) 

so that 

w 2 = 2X 0 (P pc/Z P E s f (13) 

where 

E s = V(4n p e 2 /a) = 21.2 MeV (14) 

with a the fine structure constant, p the ion momentum 
array, the array fi of ion speeds relative to the speed of light 
c, Z P the array of projectile charges, and X 0 is the electron 
radiation length in the material. The electron radiation 
length (g/cm 2 ) is given by 

Xq 1 = 4a (N a /A) Z T 2 r 0 2 ln(181Zj m ) (15) 

with Avagadro’s number N A , A the molecular weight, and r 0 
the classical electron radius. We will be using Schimmerling and coworkers modifications to Fermi’s width formula 
(1986, Wong et al. 1989). Multiple Coulomb scattering played a critical role in prior experimental validation of the 
transport solutions (Shavers et al. 1990, 1993). An example of multiple Coulomb scattering is given in Fig. 3 
showing the emerging ion angular distribution on the beam axis and off the beam axis for 600 A MeV iron ions in 
aluminum. The iron ions on the beam axis remain highly peaked in the forward direction while those that have 
scattered off axis exhibit wider angular divergence. In all cases the angular dispersion is small and clearly will be 
important only for the uncollided beam ions. Recent advances in coupling multiple Coulomb elastic scattering with 
first order transport processes (energy loss and straggling) is discussed in Mertens et al. (2007). 



Fig. 3 Multiple Coulomb scattering of 600 A MeV 
iron ions in a 5 g/cm 2 aluminum target. 



Fig. 4 Isotropic and forward neutron spectra produced 
by 500 MeV proton in aluminum. 
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The third order physical perturbation involves the nuclear reactive processes represented by the operator I r 
of equation (6), rewritten as follows 

[D-I at -I el + oJ»*={fa(Q,Q',E,E')dQ'dE}»* (16) 

The off-diagonal nuclear reactive differential cross sections can be written in the following form 

o jKr (Q,Q’,E,E') = a jkJ JE,E')/4ji + o jKfor (Q,Q',E,E r ) (17) 


where the first term is isotropic and associated with lower energy particles produced including target fragments and 
the second term is highly peaked in the forward direction and is associated mainly with direct quasi-elastic events 
and projectile fragmentation products (Wilson 1977, Wilson et al. 1988). Surprisingly, even nucleon-induced 
reactions follow this simple form and the isotropic 
term extends to relatively high energies (see Fig. 4). 

For nucleon induced reactions, the following form has 
been used in versions of FLUKA (Ranft 1980) as 
follows 


o lk ,(Q,Q',E,E') = v jk (E') o jkr (E')f jk (E,E') 
x g R (0*Q' r E,A T ) 


(18) 


where v jk ( E) is multiplicity and the Ranft factor used 
in FLUKA is 

g R (Q‘Q',E,A T ) = N r exp[-tf/X R ] Jjy2>9>0 (19) 

and constant for larger values of production angle 6 , 
N r is normalization constant, and X R given by Ranft 
as 


= (0.12 + 0.00036 A t /E) 


( 20 ) 



Fig. 5 Probability of nuclear reaction as a function of ion type 
and energy. 


although new generalized fits are being derived. This 
separation in phase space will be further exploited in 

computational procedures. The heavy ion projectile fragment cross sections are further represented by 
or k ,ja,Q',E,E')= OjJE') N, exp [-2m vf E ")(1-Q •Q')l e, jk ] x exp[ - (E + A jk -E') 2 /2 e jk 2 ]/V(2jze jk 2 ) 


( 21 ) 


where X jk is related to the momentum downshift, s jk is related to the longitudinal momentum width, s tJk is related to 
the transverse momentum width, and N t is the transverse normalizing factor. Since the transverse width is small 
compared to the projectile and fragment energy the transverse function is highly peaked about the forward direction 

(Q 0 Q'~ 1). 

Atomic interactions limit the contributions of charged particles in the transport process. For example, the 
protons and alpha particles produced in aluminum below 100 A MeV contribute to the fluence only within a few 
centimeters of their collision source and the heavier ions are even more restricted (see Fig. 1). This is an important 
factor in that the transported secondary charged particle flux tends to be small at low energies and the role of 
additional nuclear reactions are likewise limited (see Fig. 5). 

The reaction cross section is separated by equation (17) into isotropic and forward component for which 
equation (16) may be written as coupled equations 


ID - l at -l el + Orl • 9 for = { I a r/ jQ,Q',E,E') dQ'dE '}• 9 for = S r f 0r • <*> 


and 


w for 


[D - I a ,-I d + a r ] • 9 iso = {\ a(Q,Q',E,E') dQ’dE'}’ <P /( „+ {\o.JQ,Q',E,E) dQ'dE'}’ 9 for *B r ’ <P 


I 'iso T Sr, iso w Wfor 


( 22 ) 

(23) 


Equation (22) can be written as a Volterra equation (Wilson 1977, Wilson et al. 1991) and solved as 


9 for - [G + G •S r ,for m G + G m S r for # G m 3 r ,for # G +. . ./ •&[} (24) 

for which the series can be evaluated directly as described elsewhere (Wilson et al. 1994a) and is amenable to 
numerical marching procedures. 
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The cross term in equation (23) gives rise to an 
isotropic source of light ions of only modest energies and 
neutrons. The high-energy portion of the isotropic spectra 
arises from multiple scattering effects and the Fermi motion of 
the struck nucleons within the nucleus. The low-energy 
isotropic spectra arise from decay processes of the struck 
nucleus. Spectral contributions to the Neumann series depend 
on the particle range and probability of surviving nuclear 
reactions that establish the functional form of the G matrix. 

The second term of the Neumann series is proportional to the 
probability of nuclear reaction that is limited by the particle 
range as discussed above and shown in Fig. 5. It is clear from 
Fig. 5 that those nuclear reactions for the charged particles 
below a few hundred A MeV are infrequent for which fast 
convergence of the Neumann series results. For the moment 
we will neglect the straggling and multiple-elastic processes to 
simplify the present explanation (these provide only minor corrections to space radiation exposures (Tweed et al. 
2006) but are important in laboratory testing) and examine the remaining reactive terms of equation (22). The 
corresponding Volterra equation is given (Wilson 1977) by 

<t>j(x,Q,E) = {S/EJP/E,) <t>j(r(Q,x), Q, E r ) + ZffW 'A, P/E ')f E -X dE "dQ ' O jkJ .jQ,Q',E', E') 

X <t> k (x+[Rj(E) - Rj ( E ')]Q, Q',E ")}/ Sj ( E) Pj ( E) 

where <p/(x, Q, E) in equation (25) represents the / h component of the 
forward flux, T is the point on the boundary connected to x along -Q 
(r Q , x in Fig. 7), E y = Rj 1 [p - d + RJ, p is the projection of x onto Q, 
and d is the projection of T onto Q. Equation (24) results from the 
Neumann series solution to equation (25). In the past we have 
expanded the angular integral over O' asymptotically(Wilson and 
Khandelwal 1974) and implemented a solution as a marching 
procedure (HZETRN, Wilson and Badavi 1986), as a perturbation 
expansion (Wilson et al. 1984), and by non-perturbative approximation 
(Wilson et al. 1994a) resulting in three distinct methods to evaluate the 
first order asymptotic terms, all of which have had extensive 
experimental validation (Shavers et al. 1993, Wilson et al. 1998, Shinn 
et al. 1998). A typical lowest order asymptotic solution of equation 
(25) is shown in Figs. 8 and 9 for a high energy Fe 56 ion incident on 
PMMA. Independent of the method used to evaluate the lowest order 
term, the first correction term is found by replacing the fluence in the 
integrand of equation (25) by the lowest order asymptotic solution 
<p ko (x,Q',E") as 

<Pj(x,Q,E) = {S/EJP/EJ <Pi(r(Q,x),Q,E) + If^dE 'AjP/E % % dE "clQ ' a jkf JQ,Q\ E \E') 

X (p k Jx+[Rj(E) - R/E')]Q',Q , ,E")}/S j (E)P/E ) (26) 

where (j)j(x,Q,E) is found as an integral over the neighborhood of rays centered on Q using the lowest order 
asymptotic solution (p ko (x,Q f f E ") along an adjacent ray directed along Q'. Note that the boundary condition 
reached along -Q' enters through the lowest order asymptotic approximation and the angular integral correction in 
equation (25) is determined by the homogeneity and angular dependence of the space radiation and radius of 
curvature of the bounding material as we have shown long ago (Wilson and Khandelwal 1974, Wilson 1977). 
These are the determinant factors of the magnitude of the first order asymptotic correction that is anticipated to be 
very small for human rated systems (large radius of curvature) in space radiation that is homogeneous and isotropic 
in most applications (Wilson et al. 1991, Wilson et al. 1994b). 

In a region of small radius of curvature the specific flux components near the site of evaluation will be 
missing contributions along adjacent rays that do not compensate losses along the ray on which the solution is 
evaluated representing the losses due to leakage. (Note, an asymptotic treatment of such small angle dependent 


(25) 



Fig. 7 Geometric relations of quantities 
useful in solving eauation f25k 


Transverse Fragment Spread 
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Fig. 6 Normalized transverse components for 
Ca fragmentation. 
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Fig. 8 Iron ion spectra after 23 g/cm 2 of PMMA. 


phenomena is the only useful approach circumventing 
large discretization errors.) This computational procedure 
is only a small addition to prior code development and 
will have little impact on computational efficiency. The 
angular dependence of the integral kernel of equation (26) 
is controlled by the forward reactive cross section 
Ojkjo/OfQ^E^E 0 ) with its highly peaked structure given 
by equations (18) or (21) depending on particle type. The 
angular dependence of the forward peak of fragmenting 
Ca ions at 100 and 1,000 A MeV is shown in Fig. 6. The 
low-energy ions with limited range have transverse 
components on the order of 10 degrees reducing to a few 
degrees at high energies. It should be clear that the added 
divergence added by multiple Coulomb scattering of such 
fragments (Fig. 3) is negligible to the large angular widths 
of the fragmentation event (Fig. 6) further justifying 
equation (25) in space applications. 

Note that the low energy ions have limited range 

and will contribute little to the transported flux (see Fig. 1) or nuclear reactions (see Fig. 5). The higher energy 
ions, with their much longer pathlengths giving more important contributions, are related to only a very small angle 
of acceptance (few degrees) at the boundary. 

The form of the kernel leads directly to a Gauss- 
Hermite expansion and evaluation over the 
angle of production. Although the neutron 
Neumann series for the forward components 
converge more slowly since their contribution to 
the neutron flux is not limited by atomic 
interactions these higher energy neutrons will be 
adequately evaluated by similar procedures. 

Higher order asymptotic terms can be evaluated 
with similar iteration of equation (26) if 
required but all indications are that the first such 
correction will be small for most space radiation 
applications. This leaves the diffuse comp- 
onents of neutrons and light ions produced from 
the collision of the forward components and 
transported by equation (23) to be resolved (see 
for example, Clowdsley et al. 2000, 2001). 

The transport of the low-energy neutron 
and light-ion isotropic sources in equation (23) 
dominates the solution below about 70 A MeV 


Fe 55 on P o ly methyl Meta chrylate (PMMA) 
(fragments only) 



1000 


Fig. 9 Iron ion fragment spectra behind 23 g/cm 2 of PMMA. 


(see Fig. 4). In this region light-ion transport is completely dominated by the atomic interaction terms and only a 
very small fraction have nuclear reactions making only minor contributions to the particle fields. This is especially 
true for the target fragments that can be solved in closed form including transition effects near boundaries of 
dissimilar materials (Wilson 1977, Cucinotta et al. 1991). The transport solution for the isotropic ion source terms 
to the lowest order perturbation is given by 

ct> Us Jx,Q,E) = S/e Ey dE 'AjP j(E ')f E “f 4n dE"dO' O jkM JQ,Q',E',E’) [R/E) - Rj(E’)]Q,Q f ,E")/ Sj(E)Pj(E) (27) 


k,for v 


and will give highly accurate solutions to equation (23) since very few of the ions will have reactions (see Fig. 5), 
but could be easily corrected using the HZETRN light-ion propagator applied to the diffusive source terms. Note 
the if "integral effectively sums the ion source terms along direction Q from the boundary to x. Also, the nuclear 
survival terms Pj(E) are all near unity (see Fig. 5 showing 1 - Pj(E)). 

The neutron transport is unhindered by atomic interaction (that is, neutron interaction with the atomic 
electrons). As a result neutron flux values are dominated by elastic and reactive nuclear processes. We further 
expand equation (23) for the single neutron component as 
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[Q*V + aj Mx,Q,E ) = f o nn (Q,Q\E,E') <p n (x,Q'F) dOt dE' + [S rMo - * for + S r > * (27 J, 


(28) 


where the last term is from coupling to the solution of 
equation (27) and is expected to make only small 
contributions to the neutron flux. In earlier work, we 
have approached solution of equation (28) using 
multigroup methods assuming a two-stream 
(forward/backward) approximation (Clowdsley et al. 
2000) with limited success. A comparison of various 
methods (fb, Pn, Sn, and Monte Carlo) of direct 
solution of equation (28) is given by Heinbockel et al 
(2003) and shown in Fig. 10. The Sn method treated 
only the neutrons below 20 MeV. One problem with 
the fb and Pn methods is the approximations tend to 
fail at the higher energies since the elastic cross 
sections in HZETRN are poorly represented at higher 
energies. Since the solution above 25 MeV is 
dominated by nonelastic processes, a perturbation 
approach is taken to circumvent limitations in the high 
energy elastic cross section representations (Wilson et 
al. 1991). 

The first term for diffusive neutron transport 
uses the lowest order perturbation similar to equation 
(27) given as 



g/cm 2 of aluminum shield exposed to the February 23, 1956 
solar particle event. 


t, iso °(x,O f E) = Efo p ~ d dx’ exp [- a n (E) x’J f E "f 4ll dE”dQ’ a jkis JQ,Q',E\E”) </> kf Jx - x’0,0’,E’’) 


(29) 


where p-d is the distance along -Q from x to the boundary, is the distance from x to the source point along - Q , 
and o n (E) is the total neutron cross section. Note that equations (27) and (29) transport all particles associated with 
the collisions of the forward component 9 for found as solution to equation (24). The second collision term 
associated with the diffuse charged particle field given by equation (27) is negligible but additional source terms 
from the lo west-order diffuse neutron solution given by equation (29) provides an additional strong source of 
diffuse neutrons. The added transport of these neutrons is given by 


[Q*v + °n(E)]d J(x.O,E) = I a (Q,Q ',E,E )dQ’dE’(p J(x,Q',E') + \a (Q,Q',E,E')dQ'dE’6 t °(x,Q f ,E') (30) 


where <\> ni J( x,Q,E) is the remaining diffuse neutron component. The source term to the far right of equation (30) 


have been solved in exact 3D geometry and the energy 
spectrum is much degraded for the source term in equation 

(29) . Typical spectra of the sources in equations (29) and 

(30) are compared in Fig. 11. The spectral properties of 
the second collision source terms of equation (30) result 
from terms like 

fi(E) = Jf(E,E') f iso (E', 500 MeV) dE' (31) 

It is clear that the diffuse spectra of neutrons from this 
term are highly degraded in energy for which Sn, Pn, or fb 
approximations are fully applicable and our final attention 
turns to solution of equation (30). Note that the software 
prepared for equation (30) is also applicable to nuclear 
reactor shielding issues with appropriate source terms. 
This is the typical nuclear engineering problem for which 
a multitude of methods have been developed such as the 
Sn, multigroup, and collocation methods already applied 
to HZETRN. It is mainly a question of computational 
efficiency and we continue to investigate this issue. 



Figure. 11. Second collision neutron spectrum fl(E). 
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Although, equation (30) exhibits behavior 
similar to thermal diffusion there are strong 
differences between thermal and neutron diffusive 
processes. Thermal diffusion at ordinary temperatures 
has minor leakage through near boundaries since 
radiative processes are proportional to T 4 (in the 
absence of convection) leaving lateral diffusion an 
important process at low temperatures. In distinction, 
neutron diffusion is dominated by leakage at near 
forward and backward boundaries in penetrating 
spacecraft walls and lateral diffusion plays a minor 
role. Generally, low-energy neutron leakage is a 
dominant process within 15-20 g/cm 2 of the bounding 
surface in most materials. Since human rated systems 
have shielding of large radius of curvature and small 
thickness (4-20 g/cm 2 ) to radius ratio as determined by 
living and working space requirements, it 
approximates a connected system of flat plates for 
which leakage at forward and backward boundaries 
dominates the transport. In this limit the transport 
simplifies to a connected set of ID transport problems 
with leakage at the back and forward boundaries 
(Clowdsley et al. 2000, 2001). 

It is well established that the light ions 
produced in spacecraft materials are well described by 
current nuclear models so long as cluster processes are 
included as shown in Figs. 12 and 13. Curve B in the 
figures include the direct processes leading to light 
ions including direct knockout and de-excitation 
processes while curve A includes cluster processes 
with comparisons to Badhwar’s PHIDE detector on 
Shuttle flight STS-48. Such solutions for charged 
particles are enhanced by the transport being 
dominated by atomic interactions. Neutron transport 
includes no atomic processes and greater demands on 
knowledge of nuclear processes and transport 
descriptions are traditionally problematic (see Fig. 
10). The approach taken here of performing a first 
order Neumann perturbation as in equation (29) 
followed by a diffusion description in equation (30) is 
a reasonable next step and shows reasonable success 
in comparison with experimental flight data, see Fig. 
14 (Getselev et al. 2004). Shown in the figure are 
measurements using photo-emulsion using two 
techniques (Li-6 loading and proton recoil methods) 
and a fission foil stack interleaved with nuclear track 
detectors. Also shown in the figure are Monte Carlo 
calculations accounting only for the GCR proton 
spectra that is expected to underestimate the spectrum 
especially at energies below 100 MeV. The geometry 
of the ISS was simplified in the Monte Carlo result 
and the only guidance we had was that the shielding 
was aluminum of 20-40 g/cm 2 thickness. In our 
simulations we included the full GCR spectrum, 
trapped protons and albedo neutrons for aluminum 



Fig. 12. Observed and calculated proton spectra. 



Klnclt energy. Me V fn 

Fig. 13. Observed and calculated deuteron spectra. 

Low Energy Neutron Flux Mir and ISS averaged for shield depth 2(M0 g/cm 2 
(Averaged for 1 990 station orbit) 



Fig 14. Observed and calculated neutron spectra. 
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spheres of thickness 20, 30, and 40 g/cm 2 in the comparisons. Clearly, the comparisons in the figure are very 
encouraging but an improved geometry definition is required to assess the accuracy of current transport methods. 

FUTURE DEVELOPMENTS 

In the past there has been a bifurcation in code development with codes used in space radiation studies 
relying on marching procedures and laboratory validation depending on Green’s function solutions. The HZE 
Green’s function has been well defined and tested in laboratory validation and directly relatable to the HZE portion 
of the space radiation shielding problem but the coupling to the light-ion/neutron propagation needs further 
investigation as does the light-ion/neutron Green’s function. 

Recent studies of the HZE Green’s function have shown that straggling, energy downshifts, and energy 
dispersion are negligible in space radiation transport but have important consequences in laboratory experiments 
(Tweed et al. 2006b, Walker et al. 2006). Similarly for multiple scattering effects as discussed above and 
elsewhere (Tweed et al. 2006, Mertens et al. 2007). This allows the first order HZE Green’s function derived many 
years ago ignoring these second order effects (Wilson et al. 1991, Chun et al. 1996) to be used in developing a new 
space radiation shielding code. The first work element is to develop an adequate Green’s function based light- 
ion/neutron propagator based on the 2005 HZETRN formalism (Wilson et al. 2006) that is applicable to solution of 
equation (26) as well through substitution of only forward components in lowest order asymptotic approximation. 

We specialize to solution along a ray Q directed along the x-axis for which the Boltzmann equation within 
the continuous slowing down approximation assuming particles are produced in the forward direction (that is, 
lowest order asymptotic approximation) 

[ d x - A;‘d E S/ e ) + a] <t>.(x, E) =IJa jk (E,E') ^(x,E r ) dE' (32) 

where o. } ( E, E ') are approximated for nucleons by the multiplicities of Ranft (1980), Bertini et al. (1972), and quasi- 
elastic contribution as described by Wilson et al. (1991). An immediate problem is the near singular nature of the 
differential operator and transformation from energy to residual range coordinates as we did in developing the 
Green’s function greatly relieves this problem (Wilson et al. 1991). Unlike the Green’s function development, 
numerical procedures are simplified by introducing only a single residual range coordinate for all ions and the 
residual proton range r is used as this common coordinate as 

r =f 0 E dE'IS(E') (33) 

and residual range of other particle types is related as v 7 « r which fails at low energies corresponding to low 
residual range due to electron capture into atomic orbitals characteristic to each ion type. The corresponding 
transport equation is 

[ d x - Vj d r + OjC r) ] ipj(x,r) = Zj r °° ( vj/v k )s jk ( r>r ' ) ') dr ' (34) 

where scaled flux is now (Vj for neutral particles such as neutrons are taken as unity in scaling relations, Wilson et 
al. 1991) 

ip/x,r) = vj S(E) <pfx, E) (35) 

and the scaled differential cross sections are 

s Jk (r,r') = S(E)O jk (E,E') (36) 

We will use the Voltera equation as derived by Wilson et al. (1991) arrived at by inverting the differential operator 
(Wilson et al. 1977, 1989,1991) of equation (34) as 

ipj(x,r) = exp[-£/r,x)] ipj(0,r + Vj x) + Ej 0 x f r + vjx "exp[-^(r,x') ](v/vj s jk ( r +Vjx' >r ') ip k (x -x',r') dr'dx' (37) 

where the exponential is the integrating factor related to attenuation of the j type ions with 

C J (r,x)=/ 0 x cr j (r + Vjx')dx' (38) 

and is the particle extension coefficient. The Volterra equation (37) can be solved either as a Neumann series or by 
numerical marching procedures. Note that the inverse mapping is taken as 

<t>fx,E) - A. ipj(x,r)IS^E) (39) 
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to guarantee the equilibrium solution given as equation (36) at low energies away from the boundaries (note, the 
proton stopping power is used in case of unsealing the neutron flux). 

The Green’s function associated with equation (34) relates the solution ffx,r) to the boundary condition 
ipj(0,r) as given by 

ipj(x,r) = I k fGj Jx,r,r ') </>/' 0,r ') dr ' (42) 

resulting in the following equation for the Green’s function 

R - v, d r + Oj(r) ] GJx,r,r ') = Zj r °° (v/vjs jk(r ,r " ) G k Jx,r 'r ") dr " (43) 

where G j>k (x y r,r') approaches a delta function in r - r' and j,k at the boundary. Functionally G jk (x,r,r') satisfies a 
Boltzmann equation as does ip/x,r) and can be solved by similar procedures. 

FUTURE NEEDS 

From a transport perspective, the individual physical processes are well developed but there remain some 
tasks to smoothly integrate these methods into well-defined and efficient computational procedures. Primary 
uncertainties remain in the nuclear cross sections required to support the transport capability. An adequate neutron 
scattering database for which the above solution techniques were designed to minimize the impact of poor cross 
section definition above 25 MeV. Even below 25 MeV, the cross sections currently integrated into the transport 
procedures are of limited accuracy (Wilson et al. 1991) and major improvements can be made in both total and 
differential cross sections. The fragmentation cross sections are mainly based on the NUCFRG2 and an upgrade to 
more accurate quantum mechanical models is necessary. Continued testing against laboratory and spaceflight 
measurements is a necessary requirement for progress. 

CONCLUSIONS 

The current approach to solving the Boltzmann equation has been developing over 30 years following on 
the defining papers of Wilson and Khandelwal (1974) and Wilson (1977) with a major milestone in the 
development in the comprehensive report by Wilson et al. (1991). At this point in the development, an array of 
methodologies have been defined to treat the coupled solutions in a consistent way as to allow approach to a final 
development meeting all of the requirements of NASA in both space exposures to protect the astronauts and 
equipment on future missions as well as defining laboratory ion beam qualities supporting NASA radiobiological 
studies. 

REFERENCES 

Armstrong, T.W., Colburn, B.L. Predictions of secondary neutrons and their importance to radiation effects inside 
the international space station. Radiat. Meas. 33: 229-234; 2001. 

Badhwar, G.D., Patel, J.U., Cuccinotta, F.A., Wilson, J.W., Measurements of the secondary particle energy spectra 
in the Space Shuttle, Rad. Meas., 24: 129-138; 1995 
Badhwar, G.D., et al., In-flight radiation measurements on STS-60. Radiat. Meas. 26:17-34; 1996. 

Badhwar, G.D., et al. Radiat. Meas. 33:235-241; 2001. 

Bertini, H.W.. Guthrie, M.P., Culkowski, A.H. Nonelastic Interactions of Nucleons and Jt-Mesons with Complex 
Nuclei at Energies Below 3 GeV. ORNL TM-3148, 1972. 

Chun, S.Y., Khandelwal, G.S., Wilson, J.W., A Green’s function method for high charge and energy ion transport. 
Nucl. Sci. & Eng. 122: 267-275; 1996. 

Clowdsley, M.S., et al. A comparison of the multigroup and collocation methods for solving the low-energy 
neutron Boltzmann equation. Can. J. Phys. 78: 45-56; 2000. 

Clowdsley, M.S., Wilson, J.W., Shinn, J.L., Badavi, F.F., Heinbockel, J.H., Atwell, W., Neutron environment 
calculations for low Earth orbit. S AE 0 1 ICES2327, 200 1 . 

Cucinotta, F. A.; Katz, R.; Wilson, J. W.; Townsend, L. W.; Shinn, J. L.; and Hajnal, F., Biological effectiveness of 
high-energy protons: Target fragmentation. Radiation Research 127:130-137; 1991. 

Cucinotta, F. A., and R. R. Dubey, Alpha-Cluster Description of Excitation Energies in ^C i}^C, 3a) X at 2.1 A 
GeV, Phys. Rev., C50: 979-984; 1994. 


10 



Cucinotta, F. A.; Wilson, J. W.; Shinn, J. L.; and Tripathi, R. K: Assessment and Requirements of Nuclear 

Reaction Data Bases for GCR Transport in the Atmosphere and Structures, Adv. Space Res., 21 ( 12 ): 1753- 
1762 1998; 

Cucinotta, F.A., Wilson, J.W., Saganti, P., Hu, X. Kim, M.Y., Cleghom, T., Zeitlin, C.E., Tripathi, R.K., Isotopic 
dependence of GCR fluence behind shielding. Radiat. Meas. In press, 2007. 

Getselev, I., Rumin, S., Sobolevsky, N., Ufimtsev, M., Podzolko, M., Absorbed dose of secondary neutrons from 
galactic cosmid rays inside the international space station. Adv. Space Res. 34 : 1429-32; 2004. 

Heinbockel, J.H., Feldman, G.A., Wilson, J.W., Singleterry, R.C., Clowdsley, M.S., Solutions to the low energy 
neutron Boltzmann equation for space applications. SAE 03ICES339, 2003. 

Hugger, C.P., Nealy, J.E., Clowdsley, M.S., Wilson, J.W., Qualls, G.D., Atwell, W., Cucinotta, F.A., Golightly, 
M.J., Semones, E., Shavers, M.R., Preliminary Evaluation of an ISS Radiation Shielding Model. AIAA 2003- 
6620,2003. 

Mertens, C.J., Wilson, J.W., Walker, S.A., Tweed, J., Coupling of multiple Coulomb scattering and energy loss and 
straggling in HZETRN. Adv. Space Res. In press, 2007. 

Nealy, J.E., Cucinotta, F.A., Wilson, J.W., Badavi, F.F., Dachev, Ts.P., Tomov, B.T., Walker, S.A., De Angelis, G., 
Blattnig, S.R., Atwell, W., Pre-engineering spaceflight validation of environmental models and the 2005 
HZETRN simulation code. Adv. Space Res. In press, 2007. 

Norbury, J.W., Dick, F., Tripathi, R.K., Eta-meson production in proton-proton collisions. NASA TP-in press, 
2007. 

Pinsky, L., Carminati, F., Ferrari, A. Simulation of space Shuttle neutron measurements with FLUKA. Radiat. 
Meas. 33 : 335-339; 2001. 

Ranft, J. The FLUKA and KASPRO Hadronic Cascade Codes. Computer Techniques in Radiation Transport and 
Dosimetry, W. R. Nelson and T.M. Jenkins, eds, Plenum Press, pp. 339-371, 1980. 

Schimmerling, W., Rapkin, M., Wong, M., Howard, J. The propogation of relativistic heavy ions in multielement 
beam lines. Medical Phys. 13 : 217-228; 1986. 

Shavers, M. R., Curtis, S.B., Miller, J., Schimmerling, W. The fragmentation of 670A MeV Neon-20 as a function 
of depth in water. II. One-generation transport theory. Radiat. Res. 124: 1 17-130; 1990. 

Shavers, M. R., Frankel, K., Miller, J., Schimmerling, W., Townsend, L. W., Wilson, J. W., The Fragmentation of 
670 A MeV Neon-20 as a Function of Depth in Water. III. Analytic Multigeneration Transport Theory. 
Radiat. Res., 134 ( 1 ): 1-14; 1993. 

Shinn, J.L., Cucinotta, F.A., Simonsen, L.C., Wilson, J.W., Badavi, F.F., Badhwar, G.D., Miller, J., Zeitlin, C., 
Heilbronn, L., Tripathi, R.K., Clowdsley, M.S., Heinbockel, J.H., Xapsos, M.A. “Validation of a 
Comprehensive Space Radiation Transport Code”, IEEE Transactions on Nuclear Science , 45 : 2711-2719; 
1998. 

Tripathi, R.K., Townsend, L.W., Kahn, F. Role of intrinsic width in fragment momentum distributions in heavy ion 
collisions, Phys. Rev. C., 48(4), R1775-R1777; 1994. 

Tweed, J., Walker, S.A., Wilson, J.W., Cucinotta, F.A., Tripathi, R.K., Blattnig, S., Mertens, C.J. Computational 
Methods for the HZETRN Code. Adv. Space Res. 35:194-201; 2005. 

Tweed, J., Walker, S.A., Wilson, J.W., Tripathi, R.K., Cucinotta, F.A., Badavi, F.F., An improved Green’s function 
code for HZE transport. ICES SAE 2006-01-2147, 2006. 

Walker, S.A., Tweed, J., Wilson, J.W., Cucinotta, F.A., Tripathi, R.K., Blattnig, S., Zeitlin, C., Heilbronn, L., 
Miller, J. Validation of the HZETRN Code for Laboratory Exposures with 1 A GeV Iron Ions in Several 
Targets. Adv. Space Res. 35 : 202-207; 2005. 

Wilson, J.W., Khandelwal, G.S. Proton dose approximation in convex geometry. Nucl. Tech. 23: 298-305; 1974. 

Wilson, J.W. Analysis of the theory of high-energy ion transport. NASA TN D-8381, 1977. 

Wilson, J.W., Townsend, L.W., Bidasaria, H.B., Schimmerling, W., Wong, M. Howard, J. 20-Ne depth-dose 
relations in water. Health Phys. 46 : 1 101-1 111; 1984. 

Wilson, J. W., Badavi, F. F., Methods of Galactic Heavy Ion Transport. Radia. Res. 108: 231; 1986. 

Wilson, J. W.; Townsend, L. W.; and Badavi, F. F.: A Semiempirical Nuclear Fragmentation Model. Nucl. Inst. & 
Methods B18: 225; 1987a. 

Wilson, J. W.; Townsend, L. W.; and Badavi, F. F.: Galactic Cosmic Ray Propagation in Earth's Atmosphere. 
Radiat. Res. 109 : 173-183; 1987b. 

Wilson, J.W. S. Y. Chun, W. W. Buck, and L. W. Townsend: High Energy Nucleon Data Bases. Health Phys., 
55 : 817-819; 1988. 


11 



Wilson, J.W., Lamkin, S.L., Farhat, H. Ganapol, B.D., Townsend, L.W., A hierarchy of transport approximations 
for high energy (HZE) ions. NASA TM-41 18, 1989. 

Wilson, J. W., Townsend, L. W., Schimmerling, W., Khandelwal, G. S., Khan, F., Nealy, J. E., Cucinotta, F. A., 
Simonsen, L. C., Shinn, J. L., Norbury, J. W. Transport Methods and Interactions for Space Radiations. 
NASA RP-1257, 1991. 

Wilson, J. W., Shavers, M. R., Badavi, F. F., Miller, J., Shinn, J. L., and Costen, R. C., Nonperturbative Methods in 
HZE Propagazition. Radiat. Res. 140: 241-244; 1994a. 

Wilson, J. W., Townsend, L. W., Shinn, J. L., Badavi, F. F., Lamkin, S. L., Galactic Cosmic Ray Transport 
Methods: Past, Present, and Future. Adv. Space Res. 14: (10) 841 - (10) 852, 1994b. 

Wilson, J. W., R. K. Tripathi, F. A. Cucinotta, J. L. Shinn, F. F. Badavi, S. Y. Chun, J. W. Norbury, C. J. Zeitlin, 
L. H. Hielbronn, and J. Miller, NUCFRG2, An Evaluation of the Semiempirical Nuclear Fragmentation 
Database , NASA TP-3533, 1995. 

Wilson, J. W., Cucinotta, F. A., Tai, H., Shinn, J. L., Chun, S. Y., Tripathi, R.K., Sihver, L., Transport of light ions 
in matter. Adv. Space Res. 21(12): 1763-1771; 1998. 

Wilson, J.W., Tweed, J., Tai, H., and Tripathi, R.K. A simple model for straggling evaluation. Nucl Inst & 
Methods B 194: 389-392; 2002. 

Wilson, J.W., Tripathi, R.K., Qualls, G.D., Cucinotta, F.A., Prael, R.E., Norbury, J.W., Heinbockel, J.H., and 
Tweed, J. Space Radiation Transport Methods Development, Adv. Space Res. 34: 1319-1327; 2004. 

Wilson, J.W., et al., A procedure for benchmarking laboratory exposures with 1 A GeV iron ions. Adv. Space Res. 
35: 185-193; 2005a. 

Wilson, J.W., et al., International space station: A test bed for experimental and computational dosimetry. Adv. 
Space Res. 37: 1656-1663; 2005b. 

Wilson, J.W., Tripathi, R.K., Cucinotta, F.A., Badavi, F.F., Standardized radiation shield design method: 2005 
HZETRN. SAE ICES 2006-01-2109, 2006. 

Wilson, J.W., Nealy, J.E., Dachev, Ts.P., Tomov, B.T., Cucinotta, F.A., Badavi, F.F., De Angelis, G., Atwell, W., 
Leutke, N., Time serial analysis of the LEO environment within ISS 6A. Adv. Space Res. In Press, 2007. 

Wong, M., Schimmerling, W., Phillips, H., Ludewigt, B.A., Landis, D., Walton, J., Curtis, S.B. The multiple 
coulomb scattering of very heavy charged particles. Medical Phys. 17(2): 163-171; 1989. 

Zeitlin, C., Heilbronn, L., Miller, J., Rademacher, S.E., Borak, T., Carter, T.R., Frankel, K.A., Schimmerling, W., 
Stronach, C.E. Heavy fragment production cross sections from 1.05 GeV/mucleion 56Fe in C, Al, Cu, Pb, 
and CH2 targets. Phys. Rev. C56: 388-397; 1997. 


12 



